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Abstract 



In order to study the tunneling of electrons through an interacting, ID, 
dimerized molecule connected to leads, we consider the persistent current in a 
ring embedding this molecule. We find numerically that, for spinless fermions, 
a molecule with a gap mostly due to interactions, i.e. a Mott-Hubbard gap, 
gives rise to a larger persistent current than a molecule with the same gap, but 
due only to the dimerization. In both cases, the tunneling current decreases 
exponentially with the size of the molecule, but more slowly in the interacting 
case. Implications for molecular electronic are briefly discussed. 



PACS numbers: 73.40.Gk, 73.61.Ph, 61.16.Ch, 72.80.Le 



Typeset using REVTgX 



1 



I. INTRODUCTION 



One of the major challenges of molecular electronics is to find molecules which could 
act as wires. For this purpose, special attention has been devoted to conjugated oligomers, 
and especially to the simplest of them, polyene, the finite size equivalent of polyacety- 
lene (C p H p+ 2)- If one considers non- interacting particles, the conductance of the set lead- 
molecule-lead is, at low bias and according to Landauer's formula |]||, proportional to the 
transmission coefficient through the molecule at the Fermi level of the leads. This approach 
was already successful in reproducing experimental results on a single C 60 molecule ||, and 
it seems natural to apply it to polyene. 

In the framework of non-interacting electrons, there is a charge gap due to the dimer- 
ization of the chain and Joachim and collaborators M found that the conductance shows a 
minimum when the Fermi energy of the leads is in the middle of this gap. At this point and 
only here, the current changes linearly with the voltage. Such a property is important to 
get a good molecular wire. Thus we have to consider the magnitude of the current at this 
point. Joachim showed that this minimum behaves as t e~' y ( 2N ~ 2 ' ) , 2N being the length of 
the molecule. The coefficient 7 grows with the dimerization, and the prefactor to depends on 
the lead- molecule contact. It turns out to be possible to tune these parameters to get non 
negligible currents in long wires (10 nm). These preliminary results were then confirmed by 
further simulations on polyenes and alkanes connected to golden leads. See also Ref. ||. 

However, there are very good reasons to believe that electronic correlations play a role 
in fixing the magnitude of the charge gap in polyene ||. We give here the main ones. 
First, the lowest excited state is dip ole- allowed, and this can be recovered theoretically 
only when correlations are included J7|. Second, the observed negative spin densities [|8]U 
can be understood only for models with interaction |T0f| - Third, photo-induced absorption 
experiments |11| showing a splitting of the soliton peak in the middle of the gap would be 
puzzling without interactions [Q. Furthermore, ab-initio calculations need to include in 



some way the correlation to get the right ground state dimerization . Therefore, we wish 
to include interactions inside the molecule and treat them exactly. In comparison with the 
free case and for a given gap, we would like to know if these correlations increase or decrease 
the lowest current. 

How can we study this current? Interesting results have been derived for a Luttinger 



liquid connected to leads [14,15], but they cannot be applied to the present case when the 



chemical potential lies in the gap. There is also a general formula describing the current 
through an interacting region connected to non-interacting leads [[UJ. However this formu- 
lation only allows an analytical evaluation if the interactions can be treated perturbatively, 
which is clearly not the case in polyene, and it is not so clear how to use numerical results 
obtained on finite size clusters for the Green's functions within this formalism. 

So, we need another way to estimate the conductance. In the non-interacting case, the 
behavior of the conductance is strongly related to the exponential decrease of the amplitude 
of the wave function at the Fermi level inside the molecule. ( We consider only one channel 
for simplicity. ) This exponential decrease also controls the persistent current in a ring 
embedding the molecule. So, persistent currents and conductance are related in this case. 
We will assume that such a relationship still holds, at least qualitatively, in the interacting 
case. More precisely, we will suppose that what happens to the minimum persistent current 
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by changing the molecular dimerization or interactions informs us on what would happen to 
the tunneling current, and we will concentrate on the persistent current because we know 
how to evaluate it numerically. 

The paper is organized as follows. We introduce the model in section II. In section III, we 
analyze in more details the relationship between the persistent current and the conductance 
in the non-interacting case. Finally we present the central results of this paper about the 
persistent currents in the interacting case in section IV. 

II. THE MODEL. 

We describe the molecule by a dimerized and interacting region of length M. It is 
embedded in a non-interacting ring of length N describing the electrode (see Figurel). 
The ring is pierced by an Aharonov-Bohm flux <£>. Denoting by L = N + M the total 
number of sites, the persistent current I($o) is related to the groundstate energy E(<$>) by 
/(<E>o) = ^MU=0o- f° r simplicity we will concentrate on the mean value of |/| over $ given 
by AI = ±\E(n) - E(0)\. 

We will use the following hamiltonian : 

H=H e + H a + H m + V + e m (1) 

with 

N-2 

H e = -te i(t > ]T c\ +lCi + h.c. 

i=0 

H a = (-ae-^c j L ^c - ae^ N c N _^) + h.c. 

M/2-1 M/2-1 

H m = —tie l<t> ^2 c N+2 i+i c N+2i + h- c - ~ ^2 c N+2i c N+2i-i + h- c - 

i=0 i=l 

L-2 

V = V niUi+i 

i=N 
i=N 

In these expressions, the operator cj- creates a fermion at site j, and the density 

operator at site i. The metallic electrode is described by a tight-binding hamiltonian H e 
with a hopping integral t. 

The molecule is described by the sum of three terms: H m is a tight-binding hamiltonian 
that describes the kinetic energy inside the molecule. It involves two alternating hopping 
integrals t\ and t 2 to take the dimerization into account. V represents the nearest neighbour 
repulsion between the particles, and e m fixes the chemical potential of the molecule. 

The molecule and the electrode are connected by a transfer term H a that allows the 
electrons to hop from one to the other. 

Finally all the hopping integrals are multiplied by a phase factor e 1 ^ with = $/L to 
describe the Aharonov-Bohm flux. 



3 



Of course a realistic calculation should deal with electrons, i.e. fermions with spin. 
It turns out however that the sizes one can reach with electrons are too small to allow a 
meaningful finite size analysis (see below). So we have decided to restrict ourselves to spinless 
fermions. This should be useful as a first step toward more realistic systems because the 
physics of Mott-Hubbard insulators - i.e. systems where the charge gap is due to correlations 

— is very similar for spinless fermions and fermions with spin. 

III. RELATIONSHIP BETWEEN THE PERSISTENT CURRENT AND THE 
CONDUCTANCE FOR NON-INTERACTING PARTICLES 

There is no exact relation between the persistent current and the conductance in general, 
but we will show that the two quantities share qualitative and quantitative features in the 
non-interacting case. 

A. How to calculate the conductance G? 

According to Landauer, the conductance is proportional to the transmission coefficient 
Tp at the Fermi level. This coefficient can be calculated by considering two semi infinite 
leads connected to the molecule. Using the matching procedure, the stationary state \E > 
of energy E can be searched for as: 

\E >= \k e L > +r E \ -k e L > +A™\k m > +A™\ -k m > +t E \k e R > (2) 

where r E , A™, A™ and t E are unkonwn coefficients to be determined by continuity conditions. 
Here \k L R > and \k m > correspond to the states with energy E that one would get for 
infinite leads or for an infinite molecule, and the momenta must satisfy: E = —2t cos(/c|) = 

— \Jt\ + 1t\t2 cos^fcj 1 ). The transmission coefficient is then obtained as T(E) = \t E \ 2 . 
An analytical formula for T(E) can be derived, but is is too complicated to be written down. 
It is also possible to get T(E) for more realistic band structures. See Ref. |T7|] . A typical 
example is depicted in Figure 2. Although the method is different from that of Ref. [[|, we 
have checked that the results are indeed the same. 

B. How to calculate the persistent current AI ? 

We now consider the geometry of Figure 1. In the absence of correlations, the total 
energy E inside the ring is the sum of the energies Ei of the occupied monoparticular states 
(up to E E ). Thus the total current I is the sum of the individual currents Ii of these levels. 

The simplest way to get these energies Ei is again to use the matching procedure. The 
ring being closed, a monoparticular state \Ei > must now be searched for as: 

\Ei >= A\\k\ > +A e 2 \k e 2 > +A™\k™ > +A™\k™ > (3) 

where \kj > and [kj 1 > are again free propagating waves in the electrode and 
in the molecule. Note however that they now satisfy: Ei = —2tcos(kj — <p) = 

— \Jt\ + t\ + 2tit 2 cos(2/c™ — 2(f)), j = 1, 2, because of the Aharonov-Bohm flux. 
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The continuity conditions yield 4 linear equations for the coeficients A\ and A™ 1 . We get 
a solution E^ each time the determinant of this system vanishes. This determinant can be 
evaluated numerically, and we get the monoparticular states and energies as in the preceding 
subsection. 

An alternative way consists in using the transfer matrix formalism. If a and (5 are the 
amplitudes of the incoming and outgoing waves on the left side of the molecule, the incoming 
and outgoing amplitudes on the right side a' and j3' are given by 

(P\ = ( I/** -rt/t%\ fa\ 
WJ \-r k /t k l/t k J\PJ 

where t k and r k the transmission and reflection coefficients for an incident wave vector k 
defined in the preceding subsection. The matrix T m that enters this equation is called the 
transfer matrix of the molecule. 

A similar definition for the transfer matrix T e of the electrode holds, with r k = and 
t k = e tkN . The continuity conditions can then be written: 

r«xr m (°)=«»(°) (4) 

Note that in this approach the flux is concentrated in the boundary conditions, so that 
the energy is related to the wave-vectors by the usual relations of subsection II- A. 

Denoting by 0^ the phase of t k , the condition that the corresponding determinant van- 
ishes yields: 

\t k \ cos(<E>) = cos(0 fe + kN) (5) 

This equation can be solved numerically for k, t k being calculated as in subsection A, 
and we get again the spectrum. While the matching procedure is more convenient for that 
purpose, this equation turns out to be very useful to discuss the minimum value of the 
persistent current as a function of the chemical potential e m (see subsection II-D). 



C. General shape of G and A J. 

We are now in a position to compare the behaviors of AI and T with the Fermi level E F 
of the electrode. More precisely, the important parameter is the difference E' = E F — e m 
which vanishes when Ep sits right in the middle of the molecular gap. 

We can first control E' by varying the value of e rn for a given filling of the ring. The plot 
of AI(E') for a half-filled system is given in Figure 3. Although their shape and size are 
modified, it is important to notice that the resonances are located at the same energies as 
for T(E F ). Besides, as long as a is not too large - which is the case in Figure 3 - the number 
of resonances is equal to M for both. Some of them are washed out - both for T and for AI 
- when a is of the order of t however. 

We can also vary the filling Np of the ring (which control Ep) for a given e m , say 0, and 
we get a similar result for AI as a function of E' (See Figure 3.). The only difference is that 
some resonances are now missing since the bandwidth of the electrode is smaller than that 
of the molecule in the example of Figure 3. Of course they can all be restored by tuning the 
chemical potential of the molecule e m . 
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D. Value of the minimum current 



We now compare the minimum value of the persistent current and the minimum value 
of the transmission coefficient. They both appear for Ep = 0, if the chemical potential of 
the molecule e m is set to 0. Let us first calculate the mean current of the monoparticular 

state located at the Fermi level Alp = —\Ek F ($ = it) — Ek F (Q = 0)1- The energy at the 

7T 

Fermi level is given by: 

E F (&) = -2tcos[fc F ($)] (6) 

For a finite-size system, Ep < and \Ep\ « 1, and Eq.(|6]) gives: 

E F ($) = 2tk F ($>) - tn (7) 

We also know that |t F ($)| « 1 and thus Eq. flD yield: 

Nkp(<&) + 6 F ($) = 2nn + evr/2 - ecos$|i F ($)| (8) 

where e = ±1 and n G Z depend on N and the scattering properties of the molecule, but 
do not depent on <p since kN + is a continuous function of (f) 
We get from Eq. (0) and Eq. (g) : 

E F ($) = -2t[6 F ($) + 2nvr + evr/2 - ecos($)|t F $)|)/iV - trr] (9) 

Thus from Eq.(^) the current at the Fermi level is : 

AI F = 2t\e(\t F (n)\ + \t F (0)\) - e F (vr) + Q F (0)\L/N (10) 

When N goes to infinity, N/L tends to 1, |t F (vr)| and |t F (0)| merge and (Q f (tt) — f (O)) 
goes to zero. 
Finally : 

AI F = tt^Tp~ (11) 

It is possible to see in the same way that the sign of the current carried by the successive 
levels is alternate. The absolute value of this currents changes continuously with the energy 
level, and therefore at large L, the total current, which is the sum of the currents of the 
occupied level, is half of the last one. Thus, taking the absolute value : AI = AI F /2. 

e 2 

So, by using the Landauer formula G = — T F , we get the following relation between the 

h 

persistent current and the conductance : 




e 



The amplitude of a transmitted wave after the molecule is proportional to e~ 2lM with 
2 cosh 7 = ti/t 2 . Thus T F and G decrease like e~ 2lM while AI is proportional to e~ lM . This 
point can be easily understood if we remember that we are free, by a gauge transformation, 
to concentrate the whole phase factor $ = 7r in the middle of the molecule. The difference 
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between the zero flux case and the $ = n case in the hamiltonian is SH = 2x£ 2 (cmCm_ 1 +hc). 

2 2 1 

This can be seen as a perturbation and its first order influence on the level Ei is : 

8Ei =< Vi\5H\Vi >= 2t 2 [^i(M/2) x #*(M/2 - 1) + cc] 

At the Fermi level we have : #/(M/2) oc e _7M/2 and thus 5E F oc e~ lM . This energy shift is 
proportional to the Fermi- level current and thus to the total current. 

As a conclusion, in the non-interacting case, the persistent current gives valuable infor- 
mations on the conductance. First, the resonances have the same location. Second, the 
minimum persistent current is proportional to the square root of the minimum conductance. 
Thus, since the later decreases exponentially with the size of the molecule, the former has 
the same behavior, but with a coefficient in the exponential twice as small. 

IV. PERSISTENT CURRENT FOR AN INTERACTING MOLECULE 

We now study the minimum persistent current in the ring when the interaction term 
V is included. We have used three different methods to get the groundstate energy of the 
interacting system. The last two are mainly used to check the reliability of the first one. 

A. Exact diagonalizations. 

We haved use Lanczos algorithm to diagonalize the hamiltonian and get the groundstate 
energy E Q (M, N) for small values of M and N. For a given molecule of length M, we wish 
to get rid of the effect of the finite size N of the electrode. This was achieved by performing 
a finite size scaling on E(L) with fixed M and increasing N. 

It turns out that a meaningful description fo the molecule requires at least M — 6. 
Besides, to perform this scaling, we need as many values of N as possible, and this can be 
conveniently done while keeping the density fixed only at half-filling. Then for fermions with 
spin | and for L = 16 the dimension of the Hilbert space is of the order of 1.65 10 8 , i.e. too 
large to be handled numerically. Thus, in order to get enough points to perform a reliable 
finite size scaling, we have decided to limit our study to interacting spinless fermions. In 
that case, we could go up to L = 26, in which case the dimension of the Hilbert space is of 
the order of 10 7 . 

If we plot the average persistent current versus 1/L, is appears that a quadratic law fits 
the numerical results quite well (see Figure 4) as long as the bandwidth of the electrode is 
smaller than the molecule's one. This case corresponds to a large enough density of states 
at the Fermi-energy of the electrode. To test the reliability of the fit, we have made the 
same plot up to L = 120 for non- interacting particles (see also Figure 4). The difference 
between the two extrapolated values for the two fits (up to L = 26 and up to L = 120) is 
never more than 5%. However, this is no proof that the same fit is accurate for interacting 
particles, and we now turn to alternative methods to check the extrapolation. 
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B. Small repulsion limit. 



We will now derive the persistent current to first order in the interaction V. Although 
this limit is not relevant for the case of polyene, it will be used to test the relibiality of our 
scaling law. The interaction term in the Hamiltonian is : 

M-2 

v = v 4ci4+iCi+i 

i=0 

The energy levels \Ej > of the non-interacting case can be calculated as in Section IB. 
They are non-degenerate and form a basis for the monoparticular states. If the Ej's are 
sorted in increasing order, the non-interacting ground state is given by: 

i$°>= n 4,10 > 

j<N F 

where c E creates one particle in state \Ej >. Besides, at half-filling, Np = L/2. The 
original creation operator cj can then be written: 

N 

4 = E c k <Ej\i> 

3=1 

The first order correction to the ground state energy AE =< $°|t / |<l )0 > is thus given by: 

M-2 

AE = V Yl \<i+l\E j >\ 2 \<i\E l >\ 2 (12) 

i=0 j<N F ,l<N F 

+ Yl < Ej\i + 1 >< i + l\Ei >< Ei\i >< i\Ej >) 

j<N F ,l>N F 

AE can be easily calculated numerically for periodic ((f) = 0) or antiperiodic ((f) — tt) 
boundary conditions up to quite large systems. We can see on Figure 4 that the quadratic 
scaling law is still valid up to L = 36 for small V. 

C. DMRG. 

To find the ground state energy of a one-dimensional system, an alternative method to 
the exact diagonalisation is the Density Matrix Renormalisation Group [Tj|]. While this 
method gives very precise results for open boundary conditions, it is always less accurate 
for closed boundary conditions. 

Treating the molecule exactly and considering half of the electrode as the growing block, 
we meet strong limitations. Indeed, if M — 6, adding 1 site at each junction between the 
molecule and the leads to avoid artefacts and 2 sites to let the system grow, we have to treat 
10 sites exactly. If we keep m = 100 states to describe this growing block, the dimension d 
of the global Hilbert space is roughly : d = m x 2 10 x m = 10 7 and it is difficult to do much 
better. 
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Therefore, we have only been able to get accurate results up to L = 30. We did not go 
further, because for L = 32 and the maximum available m (actually 120), the relative error 
on the current was already 5% in the non-interacting case. It is nevertheless satisfactory to 
see that the results are again consistent with the extrapolation of the exact diagonalization 
data, as shown in Figure 4. 

D. Results. 

1. Variation of AI with e m . 

The plot of AI versus e m (Figure 5) for a half-filled ring has the same shape as in the 
non- interacting case (Figure 3). The main effect of the interaction is to shift the location e™ m 
and the value I min of the minimum current. But the 6 resonances are still present, although 
their relative distances change with respect to the non- interacting case. The current I m in 
is the one we are interested in. And the corresponding e™ m does not change with L, but 
depends on V (it is roughly equal to the opposite of V): e™ n counterbalances the mean field 
effect of V. 

2. Effect ofV on I for a given gap. 

In order to study the effect of electron-electron interaction on the minimum current, we 
have to find a meaningful way to compare it to the non-interacting case. Now, in a realistic 
situation, the rough value of the hopping integrals t\ and t 2 is usually known - it is basically 
given by the total bandwidth 2(ti + 1 2 ) - but the precise value of their ratio, which together 
with the interaction term controls the charge gap A c , is not known as accurately. The gap 
itself is known quite accurately however. So we have decided to compare models with the 
same values of 2(t x + t 2 ) and A c . In a given class of models, V is then a function of t 2 /^i- 

The main result of this paper is that, for a given bandwidth and a given charge gap A c , 
the minimum persistent current I min is larger when the interaction is partially responsible 
for the charge gap than when the charge gap is due only to the dimerization. If we plot I min 
versus V for fixed A c and bandwidth, as in Figure 6 , this effect appears clearly, and it is 
quite substantial. For example, if 2(t\ + t 2 ) = 4t and A c = 3.02t we get a current twice as 
large for V = t than in the non-interacting case. Note on Figure 6 that each curve is limited 
to the right since there is a maximum value of V consistent with each given gap. 

In fact, if we plot the current versus the repulsion for given bands, it decreases from I 
to l\ by turning on the repulsion. But this repulsion also increases strongly the gap. And 
the noninteracting system having this new gap has a current much lower than I\. This is 
why the effect of the interaction for a fixed gap favours the tunneling of electrons through 
the molecule. 

3. Influence of the molecule 's length. 

The results of the previous subsection can be extended to M=8 and M=10. As mentioned 
in the Introduction, the current decreases exponentially with the length of the molecule 
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in the non-interacting case. We find that for various parameters (repulsion, bandwidth, 
dimerisation), this exponential behavior still holds for interacting particules, as shown on 
Figure 7. 

If we choose a set of parameters and consider the equivalent, non-interacting system 
(same gap, same bandwidth), this system turns out to depend only very slightly on the size. 
Therefore, the corresponding current decreases exponentially with M. And we find that the 
coefficient of this exponential is larger than in the interacting case. (On Figure 7 it is given 
by the slope of the lines.) In other words, the current is not only larger when interactions 
are present for a given size, but the exponential decrease with respect to the length of the 
molecule is slower in the interacting case as well. This is very important since it is this 
exponential decrease that puts limitations on the use of long molecules as wires. 

V. CONCLUSIONS 

In this paper, our goal was to estimate the tunneling current through an interacting 
oligomere molecule, like polyene. The interaction turned out to play an important role. In 
this molecule, the charge gap can be reproduced by a dimerization and no interaction. This 
would correspond to a band insulator in the limit of an infinite molecule. This gap can 
also be reproduced by interactions and a smaller dimerization. Then we tend toward the 
Mott-Hubbard gap limit. 

If we fix the band parameters, the charge gap increases when we turn on the interaction, 
and our numerical simulations show of course a decrease of the current. But for a given 
bandwidth, if we reduce the dimerization while we increase the interaction, in order to keep 
a constant charge gap, then the current grows. Thus a Mott-Hubbard gap is less damaging 
than a pure dimerization gap for the current. Furthermore, in the interacting case, the 
persitent current decreases exponentially with the size of the molecule, like in the non- 
interacting equivalent case, but more slowly. So, the above mentioned difference concerning 
the current between a Mott-Hubbard gap and a pure dimerization gap increases with the 
size of the molecule. 

The next step will be to study the effect of the spin to check whether our results will hold 
for a realistic molecule with spin degrees of freedom. Clearly, they can already be applied to 
compounds where the on-site repulsion is strong compared to the other energy scales, and 
polyene is not so far from this situation, since t 1 ~ 2.5eV, U ~ 11.5eV and V ~ 2AeV. 
Whether this remains true for more general systems is left for future work. 

We acknowledge very useful discussions with C. Bruder, C. Joachim and M. Magoga. 
We are especially indebted to C. Stafford for very useful explanations concerning persistent 
currents, and E. Sorensen for his help with the DMRG. The numerical simulations have 
been performed on the C94 and C98 of the IDRIS. 
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FIGURES 




FIG. 1. Geometry of the system. Each open circle corresponds to a site and the links corre- 
spond to the various hopping integrals 
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FIG. 2. Plot of the logarithm of the transmission coefficient T versus the Fermi energy Ep of 
the incident particle, in the non-interacting case. Here M = 20, t\ = 0.7t , ti = 0.56t, a = 0.35t, 
in order to compare this plot to Joachim's result [3]. 
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FIG. 3. Comparison of the transmission coefficient T and the average persistent current AI in 
the non-interacting case for M = 6, t\ = 2.6t , ti = 2.2t, a = 0.3t. (a) T versus the Fermi energy 
Ep. The on-site energy in the molecule e m was set equal to zero, (b) Variation of AI with Ep — e m 
when Ep is fixed by the band filling and e m varies. The parameters are iV = 60 and Ep = (half 
filling), (c) Variation of AI with Ep — e m when e m is fixed and Ep varies with the filling. The 
parameters are e m = and TV = 120. 
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FIG. 4. Plot of the average persistent current AI versus 1/L. (a) Non-interacting case up to 
L = 140 (diamonds), (b) First order perturbation according to Eq. [6] up to L = 36 for V/t = 0.2 
(squares), (c) Exact diagonalizations up to L = 26 (triangles) and DMRG for 28 and 30 sites (open 
circles) for V/t = 1 and e m /t = 0.76. In all cases a quadratic fit is correct. In this example M = 6, 
h = 2.6t , t 2 = 2.2t, a = 0.3t. 
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FIG. 5. Plot of the average persistent current I versus the molecular potential in the interacting 
case for V/t = 1. This curve has the same general shape as I{e m ) in the non- interacting case. Here 
M = 6, L = 18, ti = 2.6/ , t 2 = 2.2b, a = 0.3/. 




FIG. 6. Plot of the minimum of the persistent current A/ m j n versus the repulsion V for fixed 
charge gap and bandwidth. 
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FIG. 7. Plot of the logarithm of AI versus the size M of the molecule. For each interact- 
ing system (circle) the equivalent non-interacting system (same gap and band-width) is plotted 
with crosses. In each example, we get straight lines and the curve with interaction is above the 
corresponding one without interaction and decreases more slowly. 
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